#include <iostream>
#include <iomanip>

#include <sds/sfile.h>
#include <sds/sdataset.h>
#include <sds/sdataview.h>

#include <cmath>

#include "cluster.h"

using namespace sds;
using namespace std;


int main(int argc, const char *argv[]) {
  SFile srcfile(argv[1]);
  SDataView<unsigned char> src[CH];
  
  for (int k = 0; k < CH; ++k) {
    src[k] = srcfile.getDataset(k).getDataView<unsigned char>();
  }

  SFile reffile(argv[2]);
  SDataView<unsigned char> ref;
  ref = reffile.getDataset(0).getDataView<unsigned char>();

  vector<SInteger> size = ref.getSize();

  vector<Cluster> clusters(256);
  for (int i = 0; i < size[0]; ++i) {
    for (int j = 0; j < size[1]; ++j) {
      Point point;
      for (int k = 0; k < CH; ++k) {
	point[k] = src[k][i][j];
      }
      clusters[ref[i][j]] += point;
    }
  }
  cout << setprecision(4);
  for (int k = 0; k < 256; ++k) {
    if (clusters[k].size() == 0) continue;
    cout << "#" << k << ' ' << '(' << clusters[k].size() << ')' << endl;
    for (int i = 0; i < CH; ++i) {
      cout << clusters[k].expected(i) << '\t';
      cout << sqrt(clusters[k].covariance(i, i)) << '\t' << '\t';
      for (int j = 0; j < CH; ++j) {
	cout << clusters[k].covariance(i, j) << '\t';
      }
      cout << endl;
    }
    cout << endl;
  }
  return 0;
}
